N89- 24642 


COMPUTATIONAL PROCEDURES FOR 
POSTBUCKLING OF COMPOSITE SHELLS 


0. M. Stanley and C. A. Felippa 
Lockheed Palo Alto Research Laboratory 
Palo Alto, California 


SUMMARY 

A recently developed finite-element capability for genera! nonlinear shell analysis, fea- 
turing the use of three-dimensional constitutive equations within an efficient resultant- 
oriented framework, is employed to simulate the postbuckling response of an axially 
compressed composite cylindrical panel with a circular cutout. The problem is a 
generic example of modern composite aircraft components for which postbuckling 
strength (i.e.. fail- safety) is desired in the presence of local discontinuities such as 
holes and cracked stiffeners. While the computational software does a reasonable 
job of predicting both the buckling load and the qualitative aspects of postbuckling 
(compared both with experiment and another code) there are some discrepancies due 
to (1) uncertainties in the nominal layer material properties. (2) structural sensitivity 
to initial imperfections, and (3) the neglect of dynamic and local material delamina- 
tion effects in the numerical model. Corresponding refinements are suggested for the 
realistic continuation of this type of analysis. 

§1. INTRODUCTION 

Advanced composite materials, due to their superior strength-to-weight ratios 
and stiffness tailorability, have become key ingredients in the design of modern 
aerospace vehicles. However, the complex structural response associated with such 
materials coupled with the intricacy of their fabrication creates harsh requirements 
for numerical simulation. 

Specifically, a problem that is of current interest to NASA/LaRC is the determina- 
tion of the postbuckling strength of thin laminated composite shells comprising 
the “skin” of stiffened air-transport fuselages [1]. These shells are required to 
maintain safe load-carrying capability substantially beyond the point at which 
skin buckling {i.e., wrinkling) occurs. To complicate matters, aircraft fuselages 
typically feature local discontinuites, such as fasteners, stiffeners and cutouts, 
which can induce high local stress gradients that tend to delaminate the compos- 
ite material. In the presence of buckling, these local delaminations can propagate 
throughout — and hence fail — a composite structure. 

It is NASA’s ultimate goal to be able to predict such phenomena analytically. 
To this end, they have asked us to match some experimental data obtained for a 
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representative test specimen [l].* We have completed the first phase of the global 
analysis, which is described in the following sections. It employs a recently devel- 
oped finite-element shell analysis capability, featuring the use of three-dimensional 
constitutive equations within an efficient resultant-oriented framework, to deter- 
mine the buckling load and explore the postbuckling regime of an elastically in-tact 
composite model. The more ambitious global/local analysis, which involves the 
prediction of local material dclamination and its interaction with global structural 
response to determine postbuckling strength, is still in the planning phase. 

§2. PROBLEM DESCRIPTION 
§2.1 Setup 

The focal problem mentioned above is depicted in Figure 1. Shown is a moder- 
ately deep (55.6 deg) cylindrical panel with a circular cutout. The panel has a 
14 in. square planform, a 15 in. radius of curvature and a radius-to-thickness ratio 
of 150. The hole is centrally located and measures 2 in. in diameter. The compos- 
ite shell wall consists of 16 layers of unidirectional graphite fibers in epoxy-resin. 
Each layer is .0056 in. thick (for a total of .086 in.) and the layers are arranged 
in the symmetric, quasi-isotropic** stacking sequence: {±45/90/0/0/90/ ± 45} 
degrees — repeated twice. The orthotropic-elastic material properties for each 
layer are listed in the figure. Note that these properties represent nominal values 
and will require some adjustment in the sequel (see §5). 

Surrounded by a metallic test frame, the appropriate boundary conditions for the 
cylindrical panel are (i) fully clamped on the bottom edge, (ii) clamped except for 
axial motion on the top edge and (iii) simply supported on the vertical edges. 

§2.2 Experimental Results 

The test conducted by NASA consisted of statically imposing a uniform end- 
shortening, 8, to determine the load-carrying capability of the panel beyond the 
initial buckling load. Experimental results are shown in Figure 2. 

Figure 2a represents a normalized “load versus end-shortening” curve. Note that 
buckling occurs abruptly, followed by a rapid drop in the axial load. Somewhere 


* This problem was first suggested by Dr. Norman F. Knight, Jr. at NASA Langley 
Research Center. We amusingly refer to it as “Knight’s problem” both as an ac- 
knowledgment to the originator and as a reminder of the many pitfalls obstructing 
its numerical solution. 

** The term quasi-istropic refers to the fact that the resultant constitutive matrix is es- 
sentially isotropic, due to a balanced sequence of fiber angles through the thickness. 
However, in contrast to truly isotropic materials, there is some additional coupling 
between bending and twisting deformations. 



between the top and bottom of this vertical branch, which spans a period of milli- 
seconds in the experiment, delamination occurs near the hole (Fig. 2b). Although 
the delamination gets progressively worse and eventually distorts the hole, there 
is a secondary (postbuckling) stiffening branch in the load/displacement curve. 
The test was stopped at the point on this secondary branch labeled “collapse”, at 
which point extensive delamination is evident (see |l] for details). 

§3. COMPUTATIONAL APPROACH 
§3.1 Formulation 

The formulation of the governing equations and associated solution algorithms 
used to analyze the above problem is outlined in Figure 3. A detailed description 
of this approach may be found in [2]. Briefly, it employs continuum- based (CB) 
shell elements, similar to the Ahmad element [3], but extended to the nonlinear 
regime and reduced to an economical resultant form in which element stiffness and 
force operators are pre-integrated through the thickness. 

Shell Equations 

Thickness pre-integration of the CB shell equations is achieved by augmenting the 
standard, Mindlin-type [4] hypotheses of straight normals and zero normal stress 
with two additional hypotheses, namely; small transverse-shear strains and mild 
taper. As argued in [2], these additional hypotheses do not significantly alter the 
range of applicability of the original formulation. The resulting theory, which is 
expressed in terms of stress-resultants rather than pointwise (continuum) stresses, 
is referred to as a continuum-based resultant (CBR) shell formulation. 

Note that unlike earlier efforts to pre-integrate the CB shell equations, which typi- 
cally assume a constant through-thickness variation of the surface metric {e.g., [5]), 
the present CBR formulation bypasses this assumption and hence is not restricted 
to very thin shells. 

Shell Elements 

To spatially discretize the CBR shell equations, a variety of shell finite elements 
have been implemented within the above framework. However, on the basis of the 
numerical evaluation conducted in [2], only the following two shell elements were 
considered for the present analysis: (i) the nine-node Heterosis (HET) element [6], 
and (ii) a new nine-node assumed natural strain (ANS) shell-element [7j. While 
both elements are parabolically curved (Fig. 3) and use standard isoparametric 
interpolation as a starting point, each departs from the basic recipe in order to 
properly represent inextensional bending deformation for thin shells. 


In particular, the 9 HET element selectively under-integrates all stiffness and force 
terms involving membrane strains to avoid membrane “locking” , and uses a mix- 
ture of Lagrange shape functions (for rotations) and Serendipity shape functions 
(for translations) to avoid “spurious modes” otherwise evoked by reduced integra- 
tion. 

In contrast, the 9-ANS element assumes an appropriate (inextensionally accurate) 
strain field from the outset, using a modified set of Lagrange shape functions and 
employing full numerical integration throughout. Due to the fact that the strains 
are assumed in the generally non- orthogonal isoparametric coordinate basis, an 
apparent advantage of the ANS approach is its decreased sensitivity to element 
mesh distortion. 

Nonlinear Solution Procedures 

After performing an updated- Lagrange linearization of the equilibrium equations, 
a modified Newton-Raphson (NR) version of the Riks-Crisfield (RC) arc-length 
control algorithm [8] is used to trace incrementally the load-displacement curve. 
The RC procedure was adopted as a convenient means of statically traversing the 
bifurcation points that arise in shell postbuckling analysis. However, such methods 
are not foolproof, and special attention by the analyst [e.g., in the selection of 
imperfections, step-sizes and error tolerances) is often required in the vicinity of 
closely spaced or multiple bifurcation points. 

Two kinds of update procedure are required to advance the solution of the discrete 
nonlinear shell equations from one load-step to the next: (i) a kinematic update 
that “accumulates” incremental nodal displacement — both translational and ro- 
tational components; and (ii) a constitutive update that generates new element 
stresses from the corresponding displacement field. 

Only the rotational part of the kinematic update is non-trivial.* In the present 
approach, the rotation increments are used to update orthogonal triads defined at 
each shell node. The triad-update algorithm involves no trigonometric functions, 
maintains orthogonality at each step and provides a shell-oriented coordinate sys- 
tem in which the normal rotational degree of freedom may be eliminated at all shell 
nodes (except at junctures). Furthermore, once the nodal triads and reference- 
surface coordinates have been so updated, the current element configuration is 
completely and uniquely defined (see [2] and [9] for details). 

The stress update is handled via an incrementally objective algorithm [10| that 


* Since the translational components of displacement are vectorial, translation incre- 
ments are simply added to obtain total displacements and hence update the nodal 
coordinates. 



features a midpoint-rule numerical integration of rate-type constitutive equations. 
For finite-strain analysis, the constitutive algorithm additionally involves shell 
thickness updates that account for large Poisson effects. These are computed (as 
in [ll]) by recovering the normal strain increments from the constitutive equations 
via the zero normal stress (ZNS) hypothesis. 

Summary of Computational Features 

The computational features of the above approach may be summarized as follows: 


• Applicable to Both Thin and Moderately Thick Shells 

• Rotations May Be Arbitrarily Large 

• Strains May Be Large ( Except Transverse Shears ) 

• Resultant Format Yields 


Cost Savings oc Number of Layers 


• Shell Elements Are Fairly Robust 

— No Locking 

— No Spurious Mechanisms 

— Low Sensitivity to Mesh Distortion (Especially ANS Elements) 


Finally, note that while the formulation allows for finite strains and inelastic ma- 
terial behavior, the present analysis simply employs orthotropic linear- elasticity 
within each layer and does not account for the material delamination (i.e., dam- 
age) observed in the NASA experiment. 

§3.2 Implementation (NICE Software Architecture) 

The shell-element capabilities mentioned above have been implemented in a mod- 
ular fashion to facilitate research and transferral to other finite element codes. 
Some of the specific shell-element functions that are available via independent 
FORTRAN?? subroutine calls are shown in Figure 4a. A complete description of 
the shell-element software (including listings) is given in Appendix S of reference 
121 - 

The host finite- element program used by the authors is actually a network of in- 
dependently executable programs (or processors) which are coordinated via high- 
level procedures written in a mathematically oriented command-language (Fig. 
4b). In such an environment, the shell-element software is embedded in a sin- 
gle processor and the global solution algorithms are implemented eis procedures', 
examples are given in [12]. 



The software architecture (t.c., utilities) used to construct this particular analysis 
system is known as NICE (for Network of Interactive Computational Elements, 

1 12,1 3)). Due to the flexibility provided by the NICE architecture and its suitability 
for nonlinear and coupled-field problems, it is currently being explored by NASA 
as the basis for a standard generic testbed system for Computational Structural 
Mechanics (see ll-l) and other presentations therein). One of the motivating factors 
for developing such a system is the implementational complexity associated with 
a comprehensive global/local analysis of the present composite-shell problem. 

§4. FINITE-ELEMENT MODELS 

Due to the physics of the problem, a full numerical model of the test specimen is 
required (Fig. 5). The slight anisotropy emanating from the composite material 
stacking sequence is only partly responsible for the lack of available symmetry.* 
As will become apparent, the nonlinear postbuckling response is inherently non- 
symmetric due to the presence of the hole and the participation of many diverse 
mode shapes. 

Several combinations of shell-element type and mesh density were employed during 
the course of the linear (pre-buckling), stability (buckling eigenvalue) and nonlin- 
ear (postbuckling) analyses. Figure 5 shows three representative grids, involving 
300, 1500 and 5000 degrees-of-freedom, respectively. These grids correspond to 
16, 80 and 256 nine-node elements (or alternatively to 64, 320 and 1024 four-node 
elements), respectively. Note that there is intrinsic element mesh-distortion in 
these models — both in-plane and out-of-plane — due to the focus on the hole 
and the curvature of the shell. However, the elements nearest the hole (where it 
counts most) have the most regular shapes. 

The coarsest grid (l) was used to verify the modeling procedure, the finest grid 
(3) was used exclusively to check convergence of the linear and eigen solutions, 
and the intermediate grid (2) became the workhorse for nonlinear analysis. Fur- 
thermore, as little difference was observed between the 9_HET and 9-ANS elements 
(§3.1) during the early stages of analysis, the 9_HET element (which is slightly less 
expensive) is featured in the analytical results that follow. 

Boundary conditions were imposed as described in Section 2 and illustrated in 
Figure 5. To simulate end-shortening, an axial force was applied in conjunction 
with a degree-of-freedom equivalence among all axial displacements on the loaded 
edge. This was done to avoid the use of specified displacements, which tend to 
complicate the adaptive, arc-length-based nonlinear solution algorithm. 


* 


We have confirmed this via numerical experiments with an isotropic model. 



§5. LINEAR (PRE-BUCKLING) ANALYSIS 


Results for the linear pre-buckling analysis are shown in Figure 6. The defor- 
mations due to an applied axial compressive load of 22480 lbs (1 kN) are shown 
magnified by a factor of 10 in the top half of the figure -- for grid 2. The corre- 
sponding distribution of the axial stress resultant, along the panel circumfer- 
ence at mid-span is shown in the bottom half of the figure — for both grids 2 and 
3. 

Regarding the displacement solution, convergence of the axial end-shortening, 
6, was achieved with grid 2. (Grid 3 yielded less than a 1% increase in end- 
shortening.) However, the converged end-shortening solution, S = .0316, is ap- 
proximately 15% lower than the experimental value; as deduced from the linear 
portion of the experimental load-displacement curve (Fig. 2). It is presumed that 
this over-estimation of the axial stiffness is due to the uncertainty in the nominal 
material properties, a conclusion that is reinforced in §8. Thus, to compensate 
for the mismatch, the lamina principle elastic modulus (Ei) is reduced by a cor- 
responding factor in the subsequent nonlinear analysis (§7). 

From a qualitative perspective, the solution shows substantial bending deformation 
in the vicinity of the hole (Fig. 6a). This suggests that geometric nonlinearity 
may be important even at relatively low load levels and diminishes the credibility 
of linear response and buckling-eigenvalue analyses. 

Regarding the linear stress solution, note that the compressive axial resultant, 
Nxi is distributed evenly along most of the panel circumference except for a very 
localized region near the hole. While grid 2 was adequate for the displacement 
solution, grid 3 provides much more accurate resolution of this stress concentration. 
In particular, grid 2 yielded a peak stress concentration factor (SCF) of 2.8; about 
14% less than the SCF obtained with grid 3. 

As experimental data were not available for verification of the computed stresses, 
the convergence of the grid 3 stress solution was inferred by comparison with 
a closed-form (asymptotic) solution due to C.R. Steele (private communication, 
Stanford University, 1985). While the closed-form solution pertained to a purely 
isotropic panel, the linear finite-element stress solutions for isotropic and quasi- 
isotropic panels were found to be quite similar. It is also interesting that the SCFs 
for both the isotropic model (3.1) and the quasi-isotropic model (3.25) are not 
very different from the classical SCF for a flat plate with a circular hole (3.0). 

Finally, note that by strongly biasing the mesh towards the hole, it was possible 
to obtain grid-3 accuracy with grid 2 for the local stress gradients. Such biasing, 
however, was found to be unnecessary in the subsequent, globally oriented buckling 
and postbuckling analyses. 



§6. STABILITY (BUCKLING) ANALYSIS 

Figure 7 shows the first 5 buckling eigenmodes for grid 2. These results represent 
perturbations about the linear pre-buckling solution described. The eigenvalues, 
A " 1.084, 1.IOG, 1.181, 1.432, 1.582, are the ratios of the corresponding buckling 
loads to the axial load applied in the linear pre-buckling analysis. As before, grid 
2 seems to provide adequate resolution, with grid 3 giving only a 2% reduction in 
the first two eigenvalues, and a A% reduction in the remaining three. 

The following observations are important for subsequent computational purposes: 
(i) the eigenvalues are closely spaced; (ii) the eigenmodes are vastly different in 
character; (iii) there is no single form of symmetry to be exploited computation- 
ally; (iv) the first buckling mode is symmetric and bears the most resemblance 
to the linear pre-buckling solution; (iv) the second and third modes possess skew 
symmetries; (v) the fourth and fifth modes are symmetric; the latter mode fea- 
turing practically no distortion of the circular hole; and (vi) higher modes (not 
shown) look much like those for a cylindrical panel without a hole, though the 
values remain closely spaced. 

Finally, it was found that the first (t.e., critical) buckling load is approximately 
25% lower than that of an identical cylindrical panel without a hole. Hence, while 
the influence of the hole on the buckling load is only moderate [i.e., relative to the 
stress concentration factor), its influence on the buckling modes is profound. As 
we shall see, the hole has an even stronger influence on the postbuckling response. 

§7. NONLINEAR (POSTBUCKLING) ANALYSIS 

In practice, more than just a linearly converged model and an adaptive solution 
strategy were necessary to obtain a reasonable nonlinear solution. One additional 
pre-requisite was a 15% reduction in the nominal elastic modulus, Ei (from 19.6 x 
10® to 17.1 X 10®), to match the linear branch of the load versus end-shortening 
curve (as explained in §5). Another important ingredient was the specification of 
initial imperfections. In this regard, three versions of the analysis were run: (1) 
one involving no imperfections; (2) one with an imperfection amplitude of 1% of 
the thickness applied to each of the first 4 buckling modes (see Fig. 7); and (3) 
one with an imperfection amplitude of 10% of the thickness applied to each of the 
first 4 buckling modes. A discussion of the results for these three cases follows. 

§7.1 No Imperfections 

One would think that the out-of-plane “imperfections” introduced by the linear 
pre-buckling solution (Fig. 6) would be sufficient to trigger a realistic buckling 
response. However, this was not the case. With no imperfections, the computed 
solution path resembled the experiment only up to the descending part of the 
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load-displacement curve (Fig. 8a). The computed curve then rolled back onto 
itself with the stiffening branch of the postbuckling curve practically aligned with 
the pre-buckling curve. 

To gain further insight, it is useful to look at the deformation and stress histories 
portrayed in Figure 8b. Shown is a sequence of computational “snapshots” taken 
at various points (?.e., load steps) on the nonlinear load-displacement curve.* The 
variation of the axial stress resultant along the mid-span circumference is plotted 
below each frame. 

Note that the deformation starts out (at load step. 10) much like the first linear 
buckling mode (Fig. 7), with inward dimples both fore and aft of the hole, then 
articulates through the second and third modes during the initial postbuckling 
phase (steps 15-20). This rotation of the two dimples about the hole is probably 
triggered by the bending/twisting coupling inherent in the composite stacking 
sequence. The dimples continue to rotate and broaden until, at step 40, the 
pattern begins to resemble the fifth linear buckling mode. Evidently, it is this 
“locking” into mode 5 that is responsible for the excessive secondary stiffening 
in the load-displacement curve (Fig. 8a). Clearly, mode 5 is an unrealistically 
stiff one, resembling what might occur if a ring stiffener had been placed around 
the hole. This is also evident in the axial stress distribution, where the stress 
concentration has practically been shed by step 40. 

§7.2 Large Imperfections 

To avoid the unrealistic mode-5 locking observed in the preceding analysis, a fairly 
large initial imperfection was introduced in the numerical model. This was accom- 
plished by adding an equal measure of each of the first 4 buckling modes to the 
initial geometry, such that the maximum radial (x.e., shell-normal) displacement 
in each mode was equal to 10% of the panel thickness. Thus, the magnitude of 
the combined radial imperfection approached 40% of the thickness at some points. 
It is emphasized that this rather arbitrary choice of imperfections was designed 
primarily to minimize the influence of mode 5. 

The computed nonlinear response for the imperfect panel is shown in Figure 9. 
Note that the secondary stiffening branch bears more resemblance to the exper- 
iment than did the imperfection-less analysis. Unfortunately, it is also true that 
the buckling (i.e., peak) load has dropped by about 15%, and now underestimates 
the experimental buckling load by more than 20%. A heuristic explanation is 

* Note that the pre-buckling phase of the analysis appears linear with respect to the 
axial displacement, 6. However, due to the rapid growth in radial displacements 
(not shown), the analysis is actually quite nonlinear from the outset; which explains 
the relatively large number of load-steps required on the “linear” branch of Fig. 8. 


provided by considering the deformation and stress histories. 

In Figure 9b, we see that with the 10%-h imperfections, the buckling pattern leaves 
mode I almost immediately and develops an intensified inward dimple on one side 
of the panel. Consequently, the full axial load is re-distributed to the other side 
of the panel (see Nj.. plots in Fig. 9b), accounting for the reduction in both the 
buckling load and the postbuckling stiffness (fig. 9a). 

§7.3 Small Imperfections (the ‘^Bottom Line”) 

Finally, the “best of both worlds” was obtained with a 1%-h imperfection in each 
of the first 4 modes. Compared with the previous analysis, the computed load- 
displacement curve (Fig. 10) shows both an increase in the maximum load and 
a decrease in the minimum load, thus bringing the solution more in line with 
experiment. 

The improved performance, obtained by reducing the imperfections, may again 
be related to the deformation history (Fig. 10b). Here, as in the case without 
imperfections, two inward dimples develop and proceed to rotate about the hole. 
Just after buckling, however, one dimple tends to deepen while the other dimin- 
ishes, and eventually there is a double snap-through. This accounts for the double 
dip in the load-displacement curve (Fig. 10a) and seems to explain why a lower 
minimum load is obtained with the smaller imperfection. 

Still, there are some serious descrepancies between analysis and experiment, 
namely: (i) a 1% under-estimation of the buckling load, (ii) a 25% over-estimation 
of the minimum load, (iii) a 30% under-estimation of the postbuckling end- 
shortening, and (iv) a 10% over-estimation in the postbuckling stiffness. These 
will be addressed in §9. 

§8. CORROBORATION WITH ANOTHER CODE 

To support the above results, obtained via the computational procedures described 
in Section 3, parallel analyses were performed with another finite-element com- 
puter code. For this purpose, we employed the STAGS code [15], which has been 
used for more than a decade by varous government agencies and industrial firms 
(e.g., NASA and Lockheed) to analyze difficult nonlinear shell problems. Another 
reason for using STAGS is that it features finite-element computational proce- 
dures that are substantially different from those used in the present approach, 
thus adding strength to an analytical comparison. 

For the linear (pre-buckling and buckling) analysis, excellent agreement was ob- 
tained between the STAGS and NICE-based solutions. For completeness, the 
STAGS runs were performed with two radically different shell-element types, which 
both converged to the same solution as the NICE/9.HET element, albeit at slower 
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rates and from below in stiffness. (Note that while the axial pre-buckling stiffness 
of the STAGS elements converged from below, the buckling loads converged from 
above). 

In particular, the comparison included: (i) the commonly used STAGS/410 el- 
ement — a flat quadrilateral plate element based on Kirchhoff-Love theory and 
cubic membrane/bending displacement interpolation; and (ii) the less frequently 
used STAGS/422 element — a quadrilateral composed of two Kirch hoff-based tri- 
angular plate sub-elements with cubic bending interpolation and quadratic mem- 
brane interpolation. It is believed that the relatively slow convergence “from 
below” of the STAGS/410 element is due to warping sensitivity, while that of the 
STAGS/422 element is probably due to the incompatibility between membrane 
and bending displacement fields for non-flat quadrilateral element shapes.* 

For the nonlinear comparisons with N1CE/9-HET, the STAGS/410 element was 
used exclusively. The resulting load-displacement curves (for 1%-h and 10%-h 
imperfections) are summarized in Figure 11. The dashed curves represent the 
NICE/9-HET solutions and the dotted curves represent the STAGS/410 solutions; 
both were obtained with Grid 2 (Fig. 5). 

Note that although the STAGS and NICE-based solutions use different finite- 
element types, large-rotation update procedures and nonlinear solution strategies, 
the correlation is remarkable especially during the postbuckling phase. Even 
the STAGS zero-imperfection analysis (not shown) resulted in the same exces- 
sive postbuckling stiffness as displayed in Figure 8a. One other point: While 
the STAGS/410 element consistently shows about a 5% higher buckling load 
than the NICE/9-HET element, thus coming closer to the experimental peak; 
STAGS/410 is actually less accurate - with respect to discretization errors 
— than NICE/9-HET. This follows from the fact that both STAGS/410 and 
NICE/9-HET converge from above in the buckling load. This was confirmed by 
running the N1CE/9-HET element with a coarser grid, for which it too showed a 
5% higher peak. That non-converged solutions compare better with experiment 
than converged ones suggests that spatial discretization is not the only source of 
error here (e.g., see §9). 


* It is interesting to note that the STAGS/422 element was used in the related study 
conducted in [l], where it yielded a 17% more flexible linear solution, and thus agreed 
better with the linear portion of the experiment. Nevertheless, it has since been 
found that the boundary conditions were not consistently applied to the element’s 
mid-side freedoms in that analysis. By correcting this implementation error, the 
17% discrepancy with the other elements has been completely eliminated. Thus, it 
appears that the “accuracy” obtained in [l] with the STAGS/422 element is due to 
compensating errors in the nominal material properties. 


§9. CONCLUSIONS 
§9.1 Summary 

The present study may be summarized as follows: 


• PURPOSE: 

— Validate continuum-based resultant (CBR) shell formulation 
— Evaluate new shell elements 

— Gain experience with composite-shell postbuckling analysis 

• RESULTS: 

— Good agreement in pre-buckling/buckling range 
— Good qualitative agreement in postbuckling range 
— Discrepancies due to: 

Material properties 

Imperfection sensitivity 

Dynamic effects 

Delamination 

The “good” agreement obtained between the present shell-element formulation 
and experiment in the pre-buckling and buckling range was possible only after 
adjusting the nominal material properties so that the linear axial stiffnesses coa- 
lesced. The material-property modification was further justified via corroboration 
with the STAGS finite-element code, which features a substantially different com- 
putational approach. 

The “best” solution for the nonlinear response was obtained by introducing small 
(1%-thickness) imperfections corresponding to each of the first four buckling 
modes. The computed load-displacement curve (Fig. 10a), which again compared 
well with STAGS (Fig. 11), still showed major discrepancies with the experi- 
ment: The discrepancy in the buckling load (which is relatively small) may be 
due to the inadequacy of adjusting only the principal layer elastic modulus, 
rather than the complete set of orthotropic material constants. The discrepancy in 
the unloading phase of postbuckling [i.e., the computed end-shortening reversal) 
is attributed to the quasi-static approximation of what is, in reality, a dynamic 
phenomenon. Finally, the discrepancy in the stiffening phase of postbuckling is 
clearly dominated by damage, i.e., delamination observed in the experiment but 
not represented in the model. 



§9.2 Recommendations 


The goal is to eliminate the discrepancies listed with a minimum of computational 
cost and complexity. To this end, the following steps are recommended: 

1) PERFORM ADDITIONAL EXPERIMENTS. First, more experimental data 
are required to verify existing computational capabilities for composite shell post- 
buckling. For example, panel imperfections should be carefully measured and 
selected strains and overall deformation patterns should be monitored at frequent 
intervals. Additionally, an isotropic panel should be tested in order to eliminate 
the material identification problem and also to provide a standard benchmark for 
shell-element evaluation. The isotropic problem would be valuable for screening 
out geometrically sensitive elements. 

2) REFINE NONLINEAR GRID. The nonlinear analysis should be repeated 
with a finer grid (e.</.. Grid 3), as convergence in the linear regime is no guarantee 
of convergence in the nonlinear regime. Moreover, a study of modal participa- 
tion both in the imperfections and in the nonlinear response may help establish 
modeling guidelines for future analysis. 

S) INCLUDE DYNAMIC EFFECTS. Dynamic effects, which are relatively 
straightforward to incorporate, should be assessed at the first opportunity. The 
analysis could be started in a quasi-static mode, switched to an explicit transient 
response algorithm during the unstable phase, and switched again to an implicit 
algorithm during the stable postbuckling phase. 

4) REPRESENT LOCAL FAILURE (DELAMINATION). To account for the 
composite delamination mechanism, appropriate failure criteria and progressive- 
failure modeling procedures need to be developed, implemented and evaluated. 
For shell-based failure criteria, accurate stress recovery is essential. Improved 
stress resolution is required both along the surface {e.g., via stress/displacement 
iterations and adaptive refinement) and through-the-thickness , as both normal and 
transverse-shear stresses can play a dominant role in delamination. Progressive 
failure may then be simulated by methods ranging from a simplified shell model 
that selectively degrades layer properties, to a full SD analysis near the hole with an 
evolving 2D/3D transition. An intermediate approach is to “split” shell elements 
along delaminating boundaries. (See Fig. 12.) The simplest approach, however, 
has obvious implementational advantages. 

5) DEVELOP EFFICIENT GLOBAL/LOCAL ALGORITHMS. Finally, there is 
a need to reduce the cost of nonlinear analysis for such problems. The cost of the 
global analysis is dominated by the large number of iterations/steps in the linear- 
to-postbuckling transition regime. Possible approaches include the reduced-basis 


technique jl6,17j, Thurston’s method [18], and improved extrapolators. For the 
combined global/local problem, where material “properties” are changing rapidly 
during postbuckling, additional features such as line-searches [19], quasi-Newton 
stiffness updates [20[ and nonlinear substructuring may greatly improve effi- 
ciency. 

We are presently acting on recommendations (2)-(3) and will report the outcome 
in a forthcoming paper. 
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Figure 1. Knight's problem, NASA test specimen. 



Figure 2. Experimental results. 
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• Continuum-based (CB) shell equations 

- 3-D Continuum equations (equilibrium/constitutive) 

- Embed shell hypotheses (straight normals, zero normal stress) 
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• Discretize via curved "isoparametric" elements 

- Selective/reduced integration (SRI) 

- Assumed natural strain (ANS) 


• Solve nonlinear matrix equations via; 

- Linearization w.r.t. currrent configuration (UL) 

- Modified NR algorithm with adaptive (RC) strategy 

- Nodal triad updates for large rotations 


• Solve rale constitutive equations via ; 

- "Midpoint rule" incremental algorithm 

- ZNS recovery of normal strains (thickness updates) 








Figure 3. Computational approach 


USER 

Ipi^cedurtI 


PROCESSORS : 


A n^CONTROL /---A 

) L|lNTERFACE)(i^^^-t^/ 


(DATA INTERFACE) 


GLOBAL DATABASE 


NICE (Network of Interactive Computational Elements) 




mONLIN) 

lELAS 

1 


y 


b. Shell Element Processor 


Figure 4 . Implementation 






















Figure 5. Finite-element models. 



a. Linear Displacement Solution 
(Magnified by 10) 



b. Linear Stress-Resultant Solution: 
versus Y (X=L/2) 


Figure 6. Linear (prebuckling) analysis 
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Figure 9. Nonlinear (postbuckling) analysis; 10%-h imperfections, 
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Figure 10. Nonlinear (postbuckling) analysis; 1%-h imperfections 
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Figure 11. 


Corroboration with the STAGS code, postbuckling analyses. 



c = CV 


SHELL ELEMENT SHELL ELEMENT 


a. Simple: Shell -Element Material -Property Degradation 



SHELL ELEMENTS SHELL ELEMENT 


b. Intermediate: Shell-Element "Splitting" 



SOLID ELEMENTS f SHELL ELEMENT 


CONSTRAINED 

INTERFACE 


c. Complex: 3D/2D Transition 


Figure 12. Composite failure-progression models. 
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